Exact Maximal Height Distribution of Fluctuating Interfaces 



O 

o 

(N 
00 

o 

B 

> 



X3 
O 

o 



> 
o 



o 
o 



o 



X 



Satya N. Majumdar ^ and Alain Comtet ^.^ 
^Laboratoire de Physique Theorique (UMR C5152 du CNRS), Universite Paul Sabatier, 31062 Toulouse Cedex. France 
"^Laboratoire de Physique Theorique et Modeles Statistiques, Universite Paris-Sud. Bat. 100. 91405 Orsay Cedex. France 
^Institut Henri Poincare, 11 rue Pierre et Marie Curie, 75005 Paris, France 

(February 2, 2008) 

We present an exact solution for the distribution P{hm,L) of the maximal height hm (measured 
with respect to the average spatial height) in the steady state of a fluctuating Edwards-Wilkinson 
interface in a one dimensional system of size L with both periodic and free boundary conditions. For 
the periodic case, we show that P{hm,L) = L~^^^ f (hmL~^^^^ for all L > where the function f{x) 
is the Airy distribution function that describes the probability density of the area under a Brownian 
excursion over a unit interval. For the free boundary case, the same scaling holds but the scaling 
function is different from that of the periodic case. Numerical simulations are in excellent agreement 
with our analytical results. Our results provide an exactly solvable case for the distribution of 
extremum of a set of strongly correlated random variables. 
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Fluctuating interfaces are amongst the most well stud- 
ied nonequilibrium systems due to their simplicity as 
well as numerous practical applications in systems such 
as growing crystals, molecular beam epitaxy, fluctuating 
steps on metals and growing bacterial colonies [1] . While 
the past studies mostly focused on the scaling proper- 
ties of the surface roughness characterized by the average 
width of the surface height [1] , the more recent theoreti- 
cal and experimental studies have dealt with a variety of 
other important characteristics of a fluctuating interface. 
These include the distribution of the width of heights in 
the steady state [2], the statistics of first-passage events 
or persistence [3,4], the density of local maxima or min- 
ima of heights [5], the distribution of the spatially av- 
eraged height [6] as well as the distribution of height at 
any fixed point in space [7] in growing one dimensional 
Kardar-Parisi-Zhang (KPZ) interfaces [8], the cycling ef- 
fects [9], the distribution of extremal Fourier intensities 
[10] etc. 

Recently Raychaudhuri et. al. [11] studied a different 
characteristic, namely the global maximal relative height 
(MRH) (measured with respect to the spatially averaged 
growing height) of a fiuctuating interface. This is an im- 
portant observable for two principal reasons. First, it 
has important technogical significance such as in batter- 
ies where a short circuit occurs when the highest point 
of a metal surface on one electrode reaches the opposite 
one [11]. Secondly, the maximal height is an extreme 
observable measuring a rare event. While the extreme 
value statistics is well understood for a set of independent 
random variables [12], only recently physicists have been 
paying attention to the distribution of the extremum of a 
set of correlated random variables, as this question is ap- 
pearing increasingly frequently in a number of problems 
ranging from disordered systems [13] to various problems 
in computer science such as growing search trees [14] and 
networks [15]. In a fluctuating interface, the heights are 
strongly correlated and hence a knowledge of the distri- 



bution of their maximum (or minimum) would provide 
important insights into this important general class of 
extreme value problems where the random variables are 
correlated. 

In Ref. [11], the authors argued quite generally that 
the MRH hm of an interface in its stationary state in 
a finite system of size L scales as the roughness of the 
surface, hm ~ for large L, where a is the roughness 
exponent. This indicates that the normalized probabil- 
ity density (pdf) of hm has a scaling form, P{hm,L) ^ 
L"" f {hm/ L")- This was demonstrated numerically in 
[11] for a one dimensional lattice model belonging to the 
Edwards- Wilkinson (EW) universality class [16], where 
a = 1/2. Further, it was argued that the scaling func- 
tion f{x) is sensitive to the boundary conditions [11]. 

In this Letter, using simple path integral techniques 
we present an exact solution of the scaling function f{x) 
for the one dimensional EW model, both for the peri- 
odic and the free boundary conditions. For the periodic 
boundary case, we show that the scaling function f{x) is 
the so called Airy distribution function (not to be con- 
fused with the Airy function) which is the pdf of the area 
under a Brownian excursion over a unit interval and has 
been well studied in the mathematics literature [17-20]. 
We also calculate exactly the corresponding scaling func- 
tion for the free boundary condition and show that it is 
different from the periodic case. All the moments of hm 
are also computed exactly for both the boundary condi- 
tions. These results are in excellent agreement with the 
simulation results obtained by the numerical integration 
of the discretized 1-d EW equation. Our results thus pro- 
vide an exactly solvable case for the distribution of the 
extremum of a set of strongly correlated random variables. 

Our starting point is the one dimensional EW model 
[16] which prescribes a linear evolution equation for the 
height H{x, t), 
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where r]{x,t) is a Gaussian white noise with zero mean 
and a correlator, {r]{x,t)r]{x't')) = 2S{x — x')5{t — t'). 
The equation (1) has a soft (zero wave vector) mode since 
the spatially averaged height H(x,t) = H{x,t)dx/L 
keeps on growing with time (typically as \ftJV) even in 
a finite system of size L. Hence, it is useful to subtract 
this zero mode from the height and define the relative 
height, h{x, t) = H{x, t)—H{x, t) whose distribution then 
reaches a stationary state; in the long time limit in a finite 
system. Note that, by definition, 



Jo 



h{x, t)dx = 0. 



(2) 



We will see later that this constraint of zero total area 
under the relative height h plays an important role in de- 
termining the MRH distribution. All the other nonzero 
modes of h evolve identically as those of the actual height 
H. 

We first consider the periodic boundary condition, 
ft,(0) = h{L). In this case, one can decompose the 
relative height h{x,t) into a Fourier series, h{x,t) = 
J2^=-c^Km,t)e^'""''=/^. Substituting this in Eq. (1), 
one finds that different nonzero Fourier modes decou- 
ple from each other and one can easily calculate any 
correlation function. In particular, it is easy to see 
that the height h{x,t) at any given point converges to 
a stationary Gaussian distribution as t — > oo, Pst{h) = 
g-ftV2»"/V27rw2 where the width w{L) = = 
^jLjVl for all L. Moreover, one can also show that 
{dxhdx'h) S{x — x') — l/L'm the stationary state. The 
local slopes dxh are thus uncorrelated except for the over- 
all constraint due to the periodic boundary condition, 
dxdxh = that gives rise to the residual 1/L term. 
These informations can be collected together to write the 
joint probability distribution of the heights (multivariate 
Gaussian distribution) in the stationary state as, 
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(3) 



where C(L) is a normalization constant and the two delta 
functions take care respectively of the periodic bound- 
ary condition and the zero area constraint in Eq. (2). 
The constant C{L) = \/2nL^^^ can be evaluated exactly 
by integrating Eq. (3) over all heights and setting it to 
unity [21]. One can check that if one integrates out all 
the heights in Eq. (3) except at one point, one recovers 
the single point stationary height distribution mentioned 
before. 

We next define the cumulative distribution of the 
MRH, F{hm,L) = Prob[max{/i} < hm, L]. The pdf of 



the MRH is simply the derivative, P{hrm L) = ^E^jriLI^ 
Clearly F(hm,L) is also the probability that the heights 
at all points in [0, L] are less than hm and can be formally 
written using the measure in Eq. (3) as a path integral, 

F{hm,L)=C{L) I du Vh{T)e~^Jo '^^'•^-''> 

J-oo Jh(0)=u 
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h{T)dT 



I{hm,L), 
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where I{hm, L) = W^-q 9{hm — h{T)) is an indicator func- 
tion which is 1 if all the heights are less than hm and zero 
otherwise. All the paths inside the path integral prop- 
agate from its initial value /i(0) = u to its final value 
h{L) = M, where u < hm (since by definition hm is the 
maximum). A change of variable, ^(t) = hm — h{T) and 
V = hm — M in the path integral in Eq. (4) gives, 

F{hm,L)=C{L) dv Vy{r)e--Jo'^^^^-y^ X 

Jo Jy(0)=v 
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where I{hm,L) = Y[T=o^(y (''')) '^^^'^ ^ = hmL- Note 
that hm appears only through A in the delta function, 
and hence F{hm,L) = !F{A,L). In subsequent calcula- 
tions, we will keep a general A in Eq. (5) and will finally 
use A — hmL. Next we take the Laplace transform with 
respect to A in Eq. (5) and identify the quantity in- 
side the exponential as the action corresponding to a sin- 
gle particle quantum Hamiltonian, H = —■^■§^ + V{y), 
where V{y) = Xy for y > and V{y) = oo for y < 0. 
The latter condition takes care of the indicator function. 
Using the standard bra-ket notation we get. 



J^(A,L)e-'^dA = C(L) / dv <v\e-'"'^\v > 
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where Tr is the trace. Thus our problem is reduced to 
calculating just the eigenvalues of the above Hamilto- 
nian H which has only bound states and hence discrete 
eigenvalues. Solving the Schrodinger equation, one finds 
that the wavefunction (up to a normalization constant) 
is simply tpsiy) = Ai [{2Xy/^{y - E/X)] where Ai{z) is 
the standard Airy function [22] . This wavefunction must 
vanish at y = which determines the discrete eigenval- 
ues, Ek = ak)?l^2~^l^ for fc = 1,2..., where afe's are the 
magnitude of the zeros of Ai{z) on the negative real axis. 
For example, one has a\ = 2.3381 . . ., a2 = 4.0879 . . ., 
as = 5.5205 . . . etc. Upon formally inverting the Laplace 
transform in Eq. (6) and putting A = hmL we find 
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Taking derivative with respect to hm in Eq. (7) and 
making a change of variable, A = sL~^/'^, we arrive at 
our main result, P{hm,L) = L~'^/^f[hmL~^/^) for all 
L, where the Laplace transform of f{x) is given by 



/ f{x)e-'''dx = sV2^y2e 

-'0 fc=l 



(8) 



Interestingly, the right hand side of Eq. (8) turns out 
precisely to be the Laplace transform of the pdf of the 
area binder a Brownian excursion over a unit interval [19]. 
A Brownian excursion over the interval [0, 1] is simply a 
Brownian motion pinned at zero at the two ends of the 
interval and conditioned to stay positive in between. In- 
verting the Laplace transform in Eq. (8) one obtains 
f{x), known as the Airy distribution function [19], 
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where bk = 2a^/27 and U{a, b, z) is the confluent hyper- 
geometric function [22]. In Fig. 1, we have plotted f{x) 
in Eq. (9) using the Mathematica and compared it with 
the numerical scaling function generated by collapsing 
the data for 3 different system sizes obtained by numer- 
ically integrating the space-time discretized form of Eq. 
(1). Evidently the agreement is very good. 

It is easy to obtain the small x behavior of x from Eq. 
(9), since only the k = 1 term dominates as a: — > 0. Using 
U{a, b, z) ~ z~°' for large z, we get as a; ^ 0, 
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This essential singular tail near a; — > was conjectured 
in [11] based on numerics, though the exact form was 
missing. The asymptotic behavior at large x is more 
tricky to derive [23] from Eq. (9). Even the calcula- 
tion of moments from Eq. (8) is rather nontrivial. How- 
ever, it is possible to write down an exact recursion re- 
lation for the moments [18,19] and using these results, 
we get (/i;;) = M„L"/2 where Mq = 1, Mi = ^/ttJs, 
Ah = 5/12, Ms = 150r/64v^, M4 = 221/1008 etc. 
Only the second moment {h"^) = 5L/12 was computed 
before in [11] by using a different method. Finally, one 
finds that for large n, M„ ~ (71/12^)"/^. Substituting an 
anticipated large x decay of the form, /(x) ~ exp [—ax''] 
in M„ = /(x)x"dx, we get A'/„ ~ (n/a6e)"/'' for large 
n. Comparing this with the exact large n behavior of M„ 
we get a = 6 and 6 = 2, indicating /(x) ~ exp[— 6x^] as 
a; — > 00. 



There is an alternative elegant probabilistic derivation 
of the above result which we outline briefly. It proceeds 
by establishing the equivalence. 



h{x) = B{x) 
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B{T)dT, 
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where h{x) is the stationary EW interface with peri- 
odic boundary condition. B(x) is a Brownian bridge (a 
Brownian motion such that B{0) = B{L) = 0) and = 
means that the left hand side (Ihs) has the same prob- 
ability distribution as the right hand side (rhs). First, 
by construction the rhs satisfies the area constraint in 
Eq. (2). Secondly, both the Ihs and rhs of Eq. (11) 
are Gaussian variables and hence to establish the equiv- 
alence in Eq. (11), it is sufficient to show that their re- 
spective two-point correlators are identical. For example, 
one finds [21] from Eq. (1) that in the stationary state, 
{h{x)h{x')) = [LV6 - L\x - x'\ + {x - x'f\/2L for all 
L. Similarly, one can calculate the two-point correlator of 
the rhs using the representation, B{t) = x{t) — tx{L)/L, 
where x{t) is ordinary Brownian motion starting at 
x(0) = and with a correlator, (a;(r)a;(r')) = min(r, r'). 
This representation guarantees that B{0) = B{L) = 0. 
We find that the two-point correlator of the rhs is exactly 
the same as {h{x)h{x')) . This establishes the equivalence 
in law in Eq. (11) rigorously. Hence, the maximum of 
h{x) will have the same distribution as the maximum of 
the rhs of Eq. (11) which, incidentally, was computed by 
Darling in the context of statistical data analysis and he 
found [17] exactly the same Laplace transform as in Eq. 
(8). 

We next consider the free boundary condition where 
the two ends of the interface are held free. In this case, 
the joint distribution of heights in the stationary state is 
given by the same formula as in Eq. (3), except without 
the delta function 6 [/i(0) — h{L)\. This changes the nor- 
malization constant to C{L) = L. However, unlike the 
simple trace in the periodic case in Eq. (6), the Laplace 
transform in the free case turns out to be more compli- 
cated [21]. Omitting details, we find the same scaling 
as in the periodic case, P{h,n,L) = L~^^'^ f (^hmL~^^'^) , 
though the scaling function has a different Laplace trans- 
form f{s) = f{x)s-^''dx, 
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where C(a) = [J^^Aiiz)dz]''/[Ai'{-a)]'' and Ai'{z) = 
dAi {z) jdz. This Laplace transform does not seem to have 
appeared before in the mathematics literature. One can 
again express the function /(a;) in terms of the conflu- 
ent hypergeometric function [21], a Mathematica plot of 
which is shown in Fig. 1 that matches well with the nu- 
merical simulations. The small x behavior can again be 
found easily and we get. 
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where C{ai) = 3.30278..., evaluated using the Math- 
ematica. Thus the function f{x) decays shghtly faster 
as a; — !■ compared to the periodic case in Eq. (10). 
We also calculated the moments exactly [21], (/iJJ^) = 
^„L"/2 where /io = 1, A^i = \/2/7r, fi2 = 17/24, 
^3 — 123V2/140\/7r etc. We found that for large n, 
l^n ~ [rt/3e]"/2 which provides the asymptotic large x 
tail of f{x), f{x) ~ exp[— 3a;2/2] that falls off less rapidly 
than the periodic case where f(x) ~ exp[— Gx^]. 
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FIG. 1. The scaling function f{x) for the MRH distribu- 
tion for both the periodic (the top 4 curves) and the free 
(the lower 4 curves) boundary conditions. In both cases, the 
numerical curves (shown by symbols) are obtained by collaps- 
ing the data from the numerical integration of the space-time 
discretized form of the EW equation (1) for 3 system sizes 
L = 256, L = 384 and L — 512. They are compared to the 
corresponding analytical scaling functions (solid lines). 

To conclude, we note that apart from the theoretical 
interests as a solvable model, many experimental systems 
are well described by the 1-d EW equation (1). Exam- 
ples include, amongst others, the high-temperature step 
fluctuations in Si(lll)-Al surfaces [4,24] and the displace- 
ments of nonmagnetic particles in dipolar chains at low 
magnetic field [25]. Besides, the displacements of beads 
in a polymer chain with harmonic interaction (the Rouse 
model [26]) also evolve via the 1-d EW equation. Thus 
our results are relevant in these systems and it would be 
interesting to see if the MRH distribution can be mea- 
sured experimentally in such systems. 

We thank Y. Shapir, S. Raychaudhuri and C. Dasgupta 
for useful correspondence and discussions. 
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